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SUMMARY 


A new two-dimensional computer code has been developed to analyze the vis- 
cous flow around unconventional airfoils at various Mach numbers and angles of at- 
tack. The Navier-Stokes equations are solved using an implicit, upwind, finite-volume 
scheme. Both laminar and turbulent flows can be computed. A new nonequilibrium 
turbulence closure model was developed for computing turbulent flows. This two- 
layer eddy viscosity model was motivated by the success of the Johnson- King model 
in separated flow regions. The influence of history effects are described by an ordi- 
nary differential equation developed from the turbulent kinetic energy equation. The 
performance of the present code has been evaluated by solving the flow around three 
airfoils using the Reynolds time-averaged Navier-Stokes equations. Excellent results 
were obtained for both attached and separated flows about the NACA 0012 airfoil, 
the RAE 2822 airfoil, and the Integrated Technology A 153W airfoil. Based on the 
comparison of the numerical solutions with the available experimental data, it is con- 
cluded that the present code in conjunction with the new nonequilibrium turbulence 
model gives excellent results. 
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CHAPTER 1. INTRODUCTION 


In the transonic regime, the viscous flow about airfoils at angle of attack is of- 
ten complex. These flows may include a few or all of the following: shock waves, 
strong shock wave/boundary layer interactions, separation bubbles and regions of 
massive separation T]. The occurrence of these adverse flow behaviors is strongly 
influenced by the airfoil geometry. In recent years, airfoil profiles have been dramat- 
ically altered to improve airfoil characteristics 2, 3,4, 5, 6 Jj. These changes include 
leading- and trailing-edge flap deflections, blunt trailing edges, laminar flow airfoils 
and supercritical-type airfoils. Geometric features promoting adverse flow patterns 
are: 1) small leading edge radii, 2) surface and surface slope discontinuities, and 
3) nonstandard camber or thickness distributions. In the present work, airfoils with 
these geometric features have been referred to as unconventional airfoils. Some typical 
unconventional airfoil shapes are depicted in Figure 1.1. 

The occurrence of laminar flow around airfoils is an exception. Boundary layers 
on aircraft wings are predominantly turbulent. With this in mind, the objective of 
the present work is to compute the turbulent flow about unconventional airfoils at 
subsonic and transonic Mach numbers. The Navier-Stokes solver developed by Cox 
[8] is used to solve the Reynolds time-averaged equations. It is well known that the 
success of a prediction method for turbulent flows depends to a large extent on the 
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NON-STANDARD CAMBER AND 
THICKNESS DISTRIBUTIONS 


Figure 1.1: Typical unconventional airfoil shapes 
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choice of the turbulence model 9;. Although progress has been made towards the 
development of a universal turbulence model for use in turbulent flow calculations, 
the goal has not been reached yet. In many of the noncompiex flows considered in 
the past, such as attached wall-bounded flows with zero or mild pressure gradients, 
calculations based on simple algebraic turbulence models have been shown to be 
within the measurement precision of experiments. The main reason for the success 
of an algebraic eddy viscosity model is that most of the basic fluid flows used as test- 
cases are nearly in a state of self-preservation or of local equilibrium in the sense 
that generation of Reynolds stress is in complete balance with its destruction 10 . 
Turbulent flows usually take up self-preserving forms if permitted by the boundary 
conditions. Near a solid surface in the so called inner layer, the velocity scale of 
turbulence depends not only on the shear stress transmitted through the layer to the 
surface, but also on the transport of Reynolds stress. In this region, the largest stress 
containing eddies have a wavelength of the order of the distance from the surface. 
Thus the turbulence length scale in the inner layer may be taken to be proportional 
to the distance from the surface. On the other hand, the outer layer is dominated by 
large eddies [101 that have considerably longer life times. As a result, history effects 
are much more important in the outer part of a turbulent boundary layer. 

It was argued by Bradshaw and Ferriss 11] that the turbulent shear stress —xihf 
is closely related to the turbulent kinetic energy and that the latter is governed by 
the turbulent kinetic energy equation. The success of the Johnson-King model [12\ is 
a testimony to the above mentioned findings of Bradshaw and Ferriss. The Johnson- 
King model performs quite well in regions of adverse pressure gradients and separation 
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but it does not always do as well in cases where simpler equilibrium (algebraic) 
eddy viscosity models produce good results. It may be possible that this model 
overemphasizes nonequilibrium effects in near equilibrium turbulent flows. 

As a result of the above observations, an effort was undertaken to develop a new 
turbulence model from a variation of the turbulent kinetic energy equation T3\ This 
two-layer model is developed directly from the available experimental observations of 
both attached and separated flows. An ordinary differential equation, derived from 
the turbulent kinetic energy equation, is used to describe the so called history effects 
of convection and diffusion of turbulent shear stress. In view of this, the model is not 
simplv an eddy viscosity model, but also contains the desirable features of a Reynolds 
stress model. Although it bears a strong resemblance to the Johnson- King model, 
the present nonequilibrium turbulence model is unique. In addition, it is much easier 
to incorporate into a Navier-Stokes solver than the Johnson-King model since it does 
not involve any iteration. 

The performance of the new turbulence model is illustrated through compar- 
isons with experimental results for turbulent flows about three different airfoils: the 
NACA 0012 airfoil, the RAE 2822 airfoil and the Integrated Technology A 153W 


airfoil. 
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CHAPTER 2. GOVERNING EQUATIONS 


Turbulence can be regarded as the time-dependent solution of the Navier-Stokes 
equations resulting from all the nonlinear processes. These nonlinear processes include 
the instabilities within the equations and the flowfield disturbances. L nfortunately. 
the association of very small length scales of the turbulence eddies along with small 
time scales prohibit the numerical solution of the complete Navier-Stokes equations 
for complex problems on present day computers. Nonetheless, useful results can be 
obtained from numerical solutions of intelligently simplified versions of the Navier- 
Stokes equations. Currently, the main thrust of engineering calculations are based 
on the Reynolds time-averaged equations il4l. These equations can be obtained b\ 
using Favre’s mass-weighted averaging as 


dp_ 

dt 


d_ 

dx 


■(pu'j) = 0 


( 2 . 1 ) 


& d 

-SP^l) + Q— + tijP ~ ( f ij - P U l' U ] ')\ = 0 


dt 


( 2 . 2 ) 


djt 

dt 



q~j + ffUjtht - u t (f i; - pu^Ujl)] = 0 


(2.3) 


where 
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Uj(x;,0 = Uj(Xj) -r-Uj/(x;,f) 
h(xi,t) = h( Xj) + hl(xi,t) 
T(x r t) = r(ij) - Tl{x^,t) 
H( Xi ,t) = H(x t ) - 
jff(lpO = p(*i) * plixi't) 


and 


p{ x i * t ) = p(x t ) - ppx^t) 
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P u i 
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ph 
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JH 
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<9r z * 3 l 3 dxfc 

,9-,' Py. 2 t 

(9xj 3 *■? dxj. 



(2.4) 


(2.5) 
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Here p is the mean density, if z is the mean velocity in the direction of x t . Ef is the 
mean total energy per unit volume, T is the mean temperature, p is the mean pressure. 
p is the coefficient of molecular viscosity, A is the coefficient of thermal conductivity 
and h is the mean enthalpy per unit mass. It has been tacitly assumed that the 
fluctuations in p and A are negligible. In practice, the viscous terms involving the 
primed fluctuations are expected to be small and are candidates for being neglected on 
the basis of order of magnitude analysis [15] . Closure of the Reynolds time-averaged 
equations requires the specification of the turbulent stresses and heat flux terms. As in 
most turbulence models, the turbulent stresses are determined from the Boussinessq 
approximation 


du i 

~P u i ,u j' = Pt A faT 



2 - 6 . d lk 

3 G dx k 


- Sj-jpk 

3 lJ 


( 2 . 6 ) 


where it is the mean turbulent' kinetic energy per unit mass; pk = ^pu^u^. 

The apparent heat flux is related to the turbulent viscosity p{. mean flow vari- 
ables, and the turbulent Prandtl number Pr t by 


—pu^thf = Cp 


JH_d f 
Prt 9x { 


( 2.1 


Experiments confirm that Pr t (the ratio of the diffusivities for the turbulent transport 
of heat and momentum) is a well behaved function across the flow :15j. Generally 
Pr^ is close to 1 and a value of 0.9 is chosen in this report. 

The above equations are nondimensionalized with respect to freestream speed of 
sound, density, temperature, coefficient of molecular viscosity and the independent 
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variables are scaled by a characteristic length (the chord length in this case). 
Since a boundary-fitted coordinate system provides a significantly improved and easier 
means of implementing the surface boundary conditions, the foregoing equations are 
transformed from a Cartesian to a generalized nonorthogonal curvilinear coordinate 
system using the transformation 

t = t 

Z = £(*i) 

7 ? = rf( x i) 

C = C(**) 

The resulting governing equations for a two-dimensional flow are written in vector 
notation as 


OQ dF dG _ Moo , OR dS_ 

dt d£ drj Re d£ dy 1 


where 


F = yr]F - XtjG 
G = -y^F-i-x^G 
R = yr/R — x yS 
S = -y^R + x^S 
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and 


Q 

F 

G 

R 

S 

E t 

T XX 

T X y 

Tyy 

yx 


yy 


{p. pit. pv. Ef} 1 

{pu.p - pir.puv. (p~ E t ) u} T 
{pi\puv. p - pi ~ . (p - E t j t } 

T 

{0, Txx* T X y . utx x ~ <?x} 

T 

{0* Tjy. Ty y . UT X y - VTyy - <Jy} 

p[e - J [ u 2 - r“); 

2 

- p t )(' 2 ^jt/c - 2 r/ x t/^ - ) 

(^ -r pf )(^yUc - -+■ £x*£ -r Ix^'p) 

^(p ~ Pt M-Si/^ - 2??yt’77 - £x u £ - 7 x^ 77 ) 


( 2 . 10 ) 


In this study, the ratio of specific heats, -y. is taken as 1.4. The coefficient of molecular 
viscosity is approximated with Sutherland's law 


a^(l -j- 198.71° /Too) 
^ “ a 2 * 198.71°/Toc 


( 2 . 11 ) 


and the Prandtl number Pt is computed from a curve fit of the data of Keenan et al. 
:16] 


Pr{f) = 0.82 - 0.293 T + 0.178T 2 + 0.027T 3 - 0.033T 4 (2.12) 


This curve is valid in the temperature range: 
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CHAPTER 3. TURBULENCE MODEL 


Development of a New Turbulence Model for Wall Bounded Shear Flows 

The turbulent energy equation for a two-dimensional incompressible mean flow 
outside the viscous sublayer is, [17] 


.dk _ .dk 
**dx dy 


du 

— uh'f — 

dy 


d_Ds 

dy 


e = 0 


(3.1) 


where 

U/ 2 4- t'/ 2 - U'f 2 

k = . e ^ v 

and 

D s - ( ^ * kv) 

\ P 

It is believed that mass-weighted averaging accounts for much of the effects of 
compressibility on turbulence for subsonic, transonic and moderately supersonic cases 
[18] . Thus equation 3.1 is assumed to hold for compressible turbulent flows when mass 
weighted variables are employed. Furthermore, if we define 


<9u; 
dx j 
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r 


— ufvL 


r / k , 


and 



(3.2) 


equation 3.1 becomes 


m dr du dD$ r 

u v — a 1 r — a i — — + o. i — - — — 

dx dy dy dy L 


(3.3) 


where it has been taken into account that is a constant [17s. The value of aq, which 
is often referred to as the turbulent structural coefficient, decreases in the presence 
of separation [19! . Along a surface where no separation occurs. a 1 is taken as 0.3 
while for surfaces with partly separated flow, aq is assigned an average \alue of 0.25 
based on the observation that — 0.2 at separation. The length scale L defined in 
equation 3.2 is the dissipation length scale which is assumed to scale with the local 
boundary lay^r thickness in the outer part of the boundary layer, i.e.. 


L = 0.096 


(3.4) 


In the outer layer ( y > 0.2<5 approximately), the advection and diffusion terms 
in equation 3.3 become important so that the history effects on turbulence cannot be 
neglected. It was observed in this study that there is a path in the outer layer along 
which the moment of the Reynolds shear stress [yr) is maximum. This location is 
found to occur near y = 0.45<5. Along this path, equation 3.3 can be written as 


drm 

umym—z v m Tm “ a\ym T m 

Ox 



+ ym 

m 



+ a\ym 


3/2 

r m 

L 


= 0 


(3.5) 


13 


It has been tacitly assumed that the path of maximum yr is parallel to the sur- 
face which is almost exact if the boundary layer approximations are used. Here and 
in what follows, the subscript “m" denotes values where the value of yr is maximum. 
In equation 3.5. Ldi/dy can be interpreted as the square root of the ratio of turbu- 
lent shear stress to density that would result if convection and diffusion effects were 
negligible. This term is then replaced by (-utvieq) 1 ^, which is assumed to be deter- 
mined from an equilibrium eddy viscosity model. Owing to its ease in computation, 
the Baldwin- Lomax model 20! is chosen here as the required equilibrium model in 
the outer layer: 


u to*e ^ ^ c V^wake^ 


(3.6) 


where K is the Clauser constant, Cep I s additional constant and 

9 — 

F wake ~ min{y m ax F FmaxX' wk ymax F ^ dl fi F max } (3.i 

The quantities ymaxp and F ma x are determined from the vorticity function 


F(y) = y j 


! _ e - j /+/26 


where 


and 


du dv 
dy dx 



yu T 

u w 


14 


The Klebanoff intermittent}' function * and are given as 

5 = : — — - — g- (3-8) 

1. - o.oi (-'kleb yiym-axp) } 

u di f ~ ( u ~ ~ v ~ * Fmax ~ ( u “ + v )min (3-9) 

The constants appearing in equations 3.6 through 3.9 are 

K = 0.0168 

Cep — 1-6 
= °- 25 

C kleb = °- 3 

For computation of turbulent flows involving the Navier-Stokes equations, the 
determination of boundary layer thickness is not straightforward. Following Stock 
and Hasse [21], the boundary layer thickness employed here is calculated from the 
Baldwin-Lomax vorticity function using 

6 = l.9i6ymax p (3.10) 

Consequently, C'jt/eh has ^ een moc hfied an( ^ assigned a value of 0.516. The diffusion 
term in equation 3.5 is modelled in a manner similar to that used by Johnson-King 
[12] and is given by 
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1/2 

/ d£s \ _ ^r 77 on'f/ r m r mai ymax ( 3 jj) 

\ dy J m a ^( 0 . — ymax ) 9m 

Here C mQ( ^ e i is a modelling constant taken to have a value of 0.3. The derivation of 
this diffusion term is given in Appendix A. Defining g = {—uivl)m • - equation 3.5 
is now rearranged to yield 


dg_ 

dx 


2um L 


~ 9 


t’m 


r 1/2 

c model T max 


ymax 


2 ximLge 2umym 2n m (O.T<5 — ymax ) </m 


(3.12) 


Equation 3.12 can be easily solved if all the terms on the right hand side except 
g and ge are determined at the previous streamwise location. The eddy viscosity in 
the outer layer is then given by 


v t 0 


1 

g2(du/dy) m 



(3.13) 


For convenience, the boundary layer approximation of the strain rate is used. 
The kinematic eddy viscosity in the inner layer is given by 


2 r\ 2 2 ; 1 zr 

1 n . — k L) y \ tc 

L i 


(3.14) 


where k is the Karman constant (0.4), Fc is a correction applied to introduce non* 
equilibrium effects in the inner layer, and D is the damping function defined as [14. 


D = 


1 — e 


</M 


(3.15) 
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where the damping length constant is given by 


A = 26 - 


UriPw p)^ “ ^ 


and 


N = 


n.8(— 1(— : 

Pe Pw 



1/2 


i>e dp 
peuj ^ dx 


In this wall region, it is still assumed that the turbulence length scale is governed 
by the local conditions, i.e.. l i = in the fully turbulent region. Since the term g e ig 
can be interpreted as a correction applied to the turbulence velocity scale in the outer 
layer, the turbulence velocity scale correction term in the inner layer, denoted by F c . 
is postulated to have the following form 


F c 


9e:g 

l-/l 


( 3 . 16 ) 


_ ur(pu>/pe) 1/ “( r u'. T _w i )[1 ~ r tt'/ r u’ i] 

1 x-Fmax 


( 3 . 17 ) 


The details of this term are described in Appendix B. 

The eddy viscosity distribution across the entire boundary layer is then 


u t = min(u to ,u ti ) 


( 3 . 18 ) 


Wake Model 


Available experimental information on the turbulent wake of a flat plate suggest 
that the wake can be divided into three regions: the near wake, intermediate region, 
and the asymptotic far wake [22. The near wake occupies a very small region down- 
stream of the trailing edge, where the logarithmic layers of the oncoming boundary 
layers are destroyed. The intermediate region is characterized by mixing between the 
two outer layers and a resulting loss of memory of the body boundary layer. This 
region is followed by the asymptotic, small-defect, far wake in which the mean flow 
and turbulence reach a state of self preservation. 

From an examination of some of the early measurements in the wakes of cylinders. 
Schlichting deduced that eddy viscosity remains constant across the wake and is 
given by 


-^- = 0.0444 (3.19) 

U e 0 


where 



— OG 

Townsend [23] has chosen the value in equation 3.19 to be 0.032 which has been 
confirmed by Rodi [24]. This modelling constant, derived on the basis of a constant 
eddy viscosity model, is not entirely realistic in view of the marked intermittency of 
the flow over a large part of the wake. Following the approach used in Ref. ,25j, but 
also including the intermittency effects, the velocity defect in the far wake can be 


written as 
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w = W 0 f (O 


(3.20) 


where f is a function of a nondimensionalized wake half-width 


f = 


16 In 2 


v t 

U e 6 


1/2 


y 

b 


Here w 0 is the velocity defect at the wake centerline and. b is defined to be the 
wake width where the velocity defect decreases to one-half of the maximum. For 


intermittencies of the form 


1 


the velocity defect can be written as 


i v 




W 0 e 


m 

m-T 



C = 4 


v t 

U e 9 


By defining a parameter 


r = 


oo -« 2 (l-^)/C 

f e ' ' 

0 


di 




the various wake properties become a function of r. 
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In the far wake, the intermittency function suggested by Townsend 2Z in boundary 
layers is 


1 

1 - 2 ( 2 //< 5) 4 


(3.24) 


Using this intermittency. r has a value of 0.96. The asymptotic eddy viscosity then 
becomes 


^ a = 0.035 U e 6 (3.25) 

The near wake eddy viscosity is calculated following the ideas set forth by Brad- 
shaw [26], Townsend [27] and C’ebeci and Meier [28]. Denoting the wake centerline 
by y u , c ii the inner eddy viscosity formula is subdivided into two parts as shown in 
Figure 3.1. When y wc i < y < j/2> following formula is used. 

. , 1/2 

um = K{-utvt) max y c 


(3.26) 
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i/.-i mav exceed that of , at large distances from the trailing edge. In order to 
prevent this undesirable property, whenever i exceeds v t Q W - the eddy viscosity is 
taken to be Uf across the wake. 

Finally, this near-wake eddy viscosity is blended smoothly with the asymptotic 
far- wake eddy viscosity i/^ a to yield the eddy viscosity in the wake. 


v tx 


u t.c 


v t 


nearwake 


u t.a , 


; -(x-x^ e _)/(206^ e ) 


[ 3.29] 


where 


v t 


nearwa 


ke 


= mm 


max ( u ti 1 ’ u ti2 )’ u t 0 ,w } 
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CHAPTER 4. NUMERICAL METHOD 


The governing equations, which are a mixed set of hyperbolic-parabolic equations 
in time, are solved using the time-dependent approach. One of the advantages of 
this approach is that separated flow regions can be computed without anv special 
treatment. In the present study, the Navier-Stokes Solver developed by Cox $ is 
used to solve the Reynolds time-averaged equations. A brief description of the scheme 
is given below. 

This code is based on an upwind-biased, finite volume scheme. The governing 


equations are written in discrete conservation law form as 


Q n ~ l - Q? ■ = J x ,A< 

x i , j ^ i , J l iJ 

Mo o 


I ^-1 2.7 ~ . d i,j + l 2 ~ 2 


V 






Re 


' ^i-1/2 ,j ^i,j+l/2 ^i.j-1/2 ^ 

^ At; 


/J 


Tl ~f~ 1 

(4.1! 


where i and j are the indices in the ( and r/ directions, respectively, and n is the 
temporal index. The numerical fluxes F and G are approximations to the true fluxes 
at the cell faces. They are computed by a flux- difference splitting analysis based on 
Roe’s approximate Riemann solver approach [291. The viscous terms R and 5 on the 
cell faces are computed with central differences. 
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First-order Scheme 

To compute the inviscid fluxes, an exact solution to an approximate Riemann 
problem is utilized. This is achieved by writing the flux difference at a cell face in a 
nonconservative form. For example, the flux difference at (i+l/2j) is written as 


AF: 


1/2 = J. , / *Vl '2 AC ? t -l/2 

i 1/2 


(4.2) 


Here .4 is the flux Jacobian OF dQ. Conservation is maintained by Roe averaging the 
i and i — 1 flow variables in the Jacobian matrix (indicated by a tilde). The metrics in 
.4 are computed with simple differences at the cell face. The fluxes at the cell faces 
can be approximated to first order with: 


i+1/2 


A+l/2 - D i+l/2*Qi+l/2 


(4.3) 


where 

D i+ 1/2 = *i+l/2 ( A »+i/ 2 ~ A i + 1 / 2 ) R i + 1/2 

Here R and R ~ * are the right and left eigenvector matrices associated with F while 
A 4- and A - are the diagonal matrices containing the positive and negative eigenvalues 
of F . This approach results in an upwind approximation to the flux at the cell face. 
The same approach is used for the calculation of G^_ \-\!2 



24 


Second-order TVD Scheme 


The first-order scheme is extended to second-order by using an approach similar 
to that of Chakravarthy and Osher 30 For example, the second-order representation 
for F is 




1/2 


1 
2 

1 - 0 


4 

1 - o 


Fi-F ,- 1 - 

D-.i - D iIl/2 A<? i+l/2j 


D 


-1 2 A ^i-l/2 " £> i-3/2 A< ^j-3/2 


(4.4) 


where 


D- = R\~R~ l 

The accuracy parameter 0 can vary from -1 to 1. The second-order flux difference 
terms are limited to provide a total variation diminishing (T\D) scheme. Limiting 
the various flux difference terms essentially reduces the amount of mass, momentum, 
and energy convected through the cell face so that nonphysical numerical oscillations 
do not occur in the solution. Numerically, the limiting has the effect of altering the 
flux representation whenever a compression or expansion wave is encountered in the 
flow. The limiter used here is the MINMOD limiter discussed in Reference L 31] . 


= $ign(x)max{0, min{\ i |,.dj/5tgn(i)]} 


The compression parameter [3 is taken to be the maximum allowable value while 


25 


maintaining the TV D property, 
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Expanding a typical flux difference term yields 
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Implicit Algorithm 

The linearized nonconservative implicit (LNI) scheme of Yee et al. [32] is used to 
enhance the stability of the integration scheme. The scheme is first-order accurate in 
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space at the n— 1 time level and second-order accurate at n. The implicit differencing 
is fully upwind: the explicit differencing is upwind biased for a,’ =# 0. Second-order 
spatial accuracy is achieved by marching the temporally first-order accurate scheme 
in time until steady-state is reached. By lagging the Jacobian matrices to time level 
n. the resulting scheme is written as 
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The parameter 0 varies the weighting on the first-order spatial terms from fully 
explicit (0) to fully implicit (1). The Jacobian matrix N arises from the linearization 
of the viscous terms. Here, only the thin-laver terms are handled implicitly. The 
equation 4.6 is solved using an approximate factorization approach. The block penta- 
di agonal algorithm is replaced with two block tridiagonal inversions, resulting in the 
following three-step procedure for each time step. 
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Step 3: Q t r Z = 


The parameter \ is an under-relaxation factor. This is needed since the residual has 
a tendency to enter limit cycles in regions with large gradients because of the discrete 
switching of the limiter. While \ - 0.9 has negligible effect on the convergence of the 
solution, it prevents the limiter from settling into stable oscillation cycles. 

Convergence to steady-state is accelerated by using a local time-stepping ap- 
proach. The local time step is 


. [(*?■''»?) mini 1 n [(4 + y v) ij] ( .i r) 

3 v 'l - (7 - l)-V/4. 2 

The value of n varies from 0, for a constant time step, to 1, for a constant local 
Courant number. A value of approximately 5 is typical for v when n is equal to L 
Values of order 10^ are used when n is equal to 1/2. For high Reynolds number 
turbulent flow calculations, a value of one-half was found to provide the best result. 

An additional feature of the numerical method includes an implicit treatment of 
boundary conditions. The appropriate number of boundary conditions to be speci- 
fied at the farfield boundary is determined based on a localized characteristic theory. 
The farfield velocity direction is estimated by adding velocity perturbations to the 
free-stream velocity. Furthermore, nonuniform free-stream initial conditions are in- 
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corporated to improve the startup convergence characteristics of the code. 

Computational Grid 

A C-grid is employed for the present airfoil calculations. To account for blunt 
trailing edges, a wake cut is placed at the mid point of the aft surface and extended 
downstream. A typical numerical grid is displayed in Figure 4.1. The resulting 
computational domain is depicted in Figure 4.2. A representative grid for the present 
calculations had a dimension of 232 x 98 grid points. All grids were obtained with the 
elliptic grid generator GRAPE 33 : . Since the GRAPE code does not patch a grid in 
the trailing edge region of a finite-thickness trailing edge, the code was modified to 
automatically grid this region. To make efficient use of the large number of grid points 
in the tj -direction at the downstream boundary, the grid generated with GRAPE was 
replaced with an algebraic grid downstream of the trailing edge. The grid is tightl\ 
packed in the tj -direction at the airfoil surface so that y was less than 0.5 at the 
first point above the surface. In most of the calculations, the farfield boundary was 
placed approximately 15 chords from the airfoil surface. 



Figure 4.1: Representative grid for a blunt trailing edge airfoil 
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Far-Field Boundary 



Figure 4.2: Computational domain for a blunt trailing edge airfoil 
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CHAPTER 5. RESULTS AND DISCUSSIONS 


In order to validate the new turbulence model developed in this study, the flow- 
fields surrounding three airfoils: the NAC A 0012 airfoil, the RAE 2822 airfoil and 
the Integrated Technology A 153W airfoil have been computed. No adjustments were 
made to the turbulence model for these different test cases. 

Since the calculations were done for a truly two-dimensional flow, an angle of 
attack correction was applied to the experimental or geometric angles of attack to 
account for wind-tunnel wall interference. As these corrections are difficult to predict 
and somewhat arbitrary, the procedure followed here was to select angles of attack 
recommended by the experimenter. The interference effects of the wind-tunnel walls 
on the free-stream Mach number were not considered although in some cases these 
corrections may have been significant. 

NAC A 0012 Airfoil 

The NACA 0012 airfoil geometry used in this study includes a blunt trailing- 
edge. The airfoil has a thickness to chord ratio of 0.12 (located at 30% of chord) 
and a leading-edge radius of 1.58% of chord. The trailing-edge thickness is 2.1 /o 
of the maximum thickness. For this case, calculations were carried out at a single 
Reynolds number of 9 x 10® and at experimental angles of attack of 1.86° and 2.86". 
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The corresponding corrected angles of attack, as recommended by Harris 34 : , were 
taken to be 1.49° and 2.26° respectively. For both the cases, the flow was forced to 
transition at 5% of chord. Results of the computations are summarized in Table 5.1 
which consists of a tabulation of computed lift and drag coefficients compared with 
the corresponding experimental values. 

Figure 5.1 shows experimental and computed pressure distributions at a Mach 
number of 0.7. with experimental and computational angles of attack at 1.86 and 
1 . 49 °. respectively. In this case there is no separation and the computed and experi- 
mental pressure distributions 34] are in close agreement with each other. Although no 
evidence of a shock wave is seen in the experimental data, a weak shock wave appears 
in the numerical solution near the chordwise location of x/c = 0.15. This difference 
may be attributed to the uncertainties in the true angle of attack and free-stream 
Mach number. The predicted lift, drag and moment coefficients are 0.233, 0.0084 and 
0.0046. respectively, while the corresponding experimental values are 0.241. 0.00/9 
and -0.005. 

A more difficult test case with a large separation corresponding to a Mach num- 
ber of 0.799 and experimental and computational angles of attack of 2.86 and 2.26 . 
respectively, is shown in Figure 5.2. The shock location predicted by the present 
nonequilibrium turbulence model agrees closely with the experimental location. More- 
over, the correct pressure plateau behind the shock wave is captured. The pressure 
distribution on the lower surface is slightly underpredicted. The predicted lift, drag 
and moment coefficients are 0.353, 0.0368 and -0.008, respectively, which compare fa- 
vorably with the experimental values of 0.39, 0.0331 and -0.017. Figure 5.3 shows the 
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Table 5.1: Experimental and Computed Force Coefficients of NACA 0012 Airfoil 




Computation 



Experiment 


Mo c 

Re x 10 -b " 

a 

O 

C d 

a 

£l 

C d i 

0.700 

9.0 

1.49 

0.233 

0.0084 

1.86 

0.241 

0.0079 

0.799 

9.0 

2.26 

0.353 

0.0368 

2.86 

0.390 

0.0331 


skin-friction coefficient on the upper surface. No experimental skin-friction data were 
available to compare with the computed results. The purpose for the presentation of 
Figure 5.3 is to exhibit the existence of a large separation region behind the shock 
wave followed by reattachment. The performance of the present turbulence model 
is compared with those of Johnson- King and Baldwin-Lomax models in Figure 5.4. 
On the lower surface, the Baldwin-Lomax model yields the best results whereas the 
Johnson- King model performs best near the shock. Overall, the present model is seen 
to provide an improvement over the Baldwin-Lomax and the Johnson-King models. 
A basic difficulty encountered while calculating this flow was the presence of large 
amplitude unsteadiness in the solution. This behavior has also been observed by- 
several other investigators [9]. 

The shock capturing capability of the TVD scheme is apparent in the pressure 
contour plot of Figure 5.5. The shock is slightly oblique to the surface of the airfoil 
because of the rapid growth of the boundary layer. The Mach number contours and 
location of the sonic line for this case is displayed in Figure 5.6. 

RAE 2822 Airfoil 

The second airfoil analyzed was the supercritical airfoil RAE 2822 which has been 
tested extensively by Cook et al. [35] . The airfoil has a maximum thickness of 12.1% 




Figure 5.1: NACA 0012: surface-pressure coefficient at M 

G-comp — 1*49°, o-exp = F86 


= 0.7. Re — 9.0 x 10 






Figure 5.4: NAC’A 0012: effects of turbulence models on surface- pressure coefficient 

at Mo c = 0.799. Re = 9.0 > 10 6 . a CO mp = 2.26°, a eX p = 2.86 c 
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of chord and a leading-edge radius of 0.827°* of chord. Calculations were performed 
for two cases with different Reynolds numbers. Mach numbers and angles of attack. 
For both sets of freestream conditions, transition was fixed at 3% of chord. Correc- 
tions were applied to the angle of attack to account for the tunnel-wall interference 
effects. Results of the computations are summarized in Table 5.2 which consists of 
a tabulation of computed lift and drag coefficients compared with the corresponding 
experimental values. 

Figure 5.7 shows experimental and computed pressure distributions at a Mach 
number of 0.73. a Reynolds number of 6.5 x 10 6 and at experimental and com- 
putational angles of attack of 3. 19" and 2.80". respectively. In this case there is no 
separation and the computed and experimental pressure distributions 36; are in close 
agreement with each other. A shock wave is predicted near x/c = 0.54 and is in ex- 
cellent agreement with the experimental result. On the lower surface, the pressure is 
slightly overpredicted in the midsection and near the trailing-edge. The effects of dif- 
ferent turbulence models on surface pressure distribution are displayed in Figure 5.8. 
The results corresponding to the Bladwin-Lomax and the Johnson-King models were 
taken from Reference [36;. It is seen that the present nonequilibrium turbulence model 
behaves more like the algebraic Baldwin- Lomax model. Larger differences between 
model predictions and experiment are indicated by the skin- friction (Figure 5.9) and 
displacement thickness (Figure 5.10) distributions. The present model seems to give 
the best overall agreement in skin-friction. Referring to Figure 5.9, it is apparent 
that the Baldwin-Lomax model predicts weak separation at the shock wave and the 
trailing-edge, whereas the Johnson-King and the present model predict no separa- 
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Table 5.2: Experimental and Computed Force Coefficients of RAE 2822 Airfoil 


Computation j Experiment 


Mo o 

Re x 10 -6 

a 

C l 


Q 

£l 



0.730 

6.5 

1 2.80 

0.820 

0.0167 

| 3.19 

0.803 

0.0168 

0.750 

6.2 

! 2.80 

0.759 

0.0246 

! 3.19 

0.743 

0.0242 


tion. With regard to the displacement thickness distributions, all the models shown 
give similar results except near the trailing-edge where the Johnson-King model gives 
a better prediction. The predicted lift and drag coefficients are 0.82. 0.0167, re- 
spectively. while the corresponding experimental values are 0.803 and 0.0168. The 
pressure contour plot for this case is displayed in Figure 5.11. Figure 5.12 shows the 
Mach number contours and the location of the sonic line. 

The next RAE case computed involves separation and corresponds to a Mach 
number of 0.75, a Reynolds number of 6.2 x 10® and experimental and computational 
angles of attack of 3.19 c and 2.80°, respectively. The surface-pressure coefficient 
distribution is shown in Figure 5.13. The effect of turbulence model on surface- 
pressure coefficient is displayed in Figure 5.14. The shock location predicted by 
the present nonequilibrium turbulence model agrees closely with the experimental 
location. However, the shock strength is slightly overpredicted and the pressure 
distribution downstream of the shock wave indicates too much pressure recovery. 
A similar trend of too much pressure recovery is also observed with the Johnson- 
King and the Baldwin-Lomax models (Figure 5.14). Figures 5.15 and 5.16 show the 
skin-friction and displacement thickness distributions. The flow immediately after 
the shock separates for a small distance and then reattaches. Both the present model 
and the Johnson-King model yield close agreement with the experimental skin- friction 




x/c 


Figure 5.7: RAE 2822: surface-pressure coefficient at A/oc, = 0.73, Re = 6.5 x 10 




3 



Figure 5.8: RAE 2822: effects of turbulence models on surface-pressure coefficient 

at Moo = 0.73, Re = 6.5 x 10®, o.comp — 2.80°, &exp = 3-19 
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Figure 5.9: RAE 2822: skin-friction coefficient at Moo = 0.73. Re = 6.5 x 10 . 

Q-comp = 2.80". Q-exp = 3.19 1 " 



x/c 


Figure 5.10: 
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RAE 2822: displacement thickness at Moo = 0.73. Re = 6.5 x 10 . 
otcomp = 2.80 , a eX p = 3.19 










x/c 

6 

Figure 5.15: RAE 2822: skin-friction coefficient at Mo o = 0.75, Re = 6.2 x 10 

0-coTTxp 2.80 , aexp = 3.19 



Figure 5.16: RAE 2822: displacement thickness at Moo 

&comp — 2.80°. a e xp ~ 3.19 


= 0.75, Re - 6.2 x 10 
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and displacement thickness data. The predicted lift and drag coefficients are 0.759 
and 0.0246. respectively, which are in excellent agreement with the corresponding 
experimental values of 0.743 and 0.0242. Solution unsteadiness was also observed in 
this case as reported previously in References 9 and [13] . 

The pressure contour plot for this case is displayed in Figure 5.17. The shock 
is again observed to be slightly oblique to the surface of the airfoil because of the 
rapid growth of the boundary layer due to separation. Figure 5.18 shows the Mach 
number contours and the location of the sonic line. The effect of grid refinement on 
the surface-pressure coefficient is shown in Figure 5.19. 

Integrated Technology Airfoil 

A family of high speed airfoils with reduced detection characteristics were de- 
veloped by the Boeing Aerospace Company and tested at NASA Langley Research 
Center [41. The development of these airfoils involved the integration of several tech- 
nologies. As a result, these are referred to as the Integrated Technology airfoils. These 
airfoils have a characteristic thin leading-edge region, a thick belly-type midsection 
and a thin trailing-edge region with reflex curvature. Computational results of flows 
around one of these airfoils. A 153VV (see Figure 5.20) will be presented here and 
compared with the experimental results [4] . This particular airfoil was designed for a 
Mach number of 0.65 and a lift coefficient range of 0.1 to 1.0 at a Reynolds number 
of three million. It has a maximum thickness of 15.8% of chord ( at 43.8% of chord), 
a leading-edge radius of 0.868% of chord and a trailing-edge thickness of 5.2% of 


maximum thickness. 








54 


The Langley 6 - by 28-inch blowdown transonic tunnel was used to test airfoil A 
153W. The airfoil model had a four inch chord and an aspect ratio of 1.50. A total 
pressure wake survey probe located 2 ./ 5 chords downstream from the trailing-edge 
was used for the calculation of Q. Surface pressures were integrated to obtain C' n . 
The flow was tripped on both the upper and the lower surfaces at 6 % of chord. The 
experimental angle of attack mentioned in Reference [4l included corrections for lift 
interference. 

Results from two flow configurations, one subcritical and the other supercritical, 
are presented here and summarized in Table 5.3. The numerical grid contained 2-1 
92 grid cells. Twenty-five of the cells in the 77 direction were packed in one-half of the 
blunt trailing-edge region. 

The first Integrated Technology airfoil calculation corresponds to a Mach number 
of 0.602. experimental and computational angles of attack of 0.82 and a Reynolds 
number of 6 x 10®. Figure 5.21 shows the comparison between the computational and 
experimental surface-pressure coefficients. Excellent agreement between the pressure 
coefficients is observed. The predicted lift, drag and pitching moment coefficients are 
0.3986. 0.00973 and -0.0737, respectively, which compare well with the corresponding 
experimental values of 0.3651. 0.0096 and -0.0619. In this case, there is no separation 
and no supersonic regions. The pressure contour plot is shown in Figure 5.22. 

The second Integrated Technology airfoil calculation corresponds to a Mach num- 
ber of 0.725, experimental and computational angles of attack of 1.65 and a Reynolds 
number of 3 x 10®. There was a large separation region on the upper surface and 
the numerical solution indicated large amplitude unsteadiness. The flow on the up- 



Figure 5.20: Integrated Technology airfoil A 153W grid 
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Table 5.3: Experimental and Computed Force Coefficients of A 153YV Airfoil 





Computation 



Experiment 
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Re x 10 _t 

* a 

c i 

C d 
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C d 

0.602 

6.0 

0.82 

0.3986 

0.0097 

: 0.82 

0.3651 

0.0096 

1 0.725 

3.0 

1.65 

0.4394 

0.0176 

1.65 

0.4378 

0.0163 


per surface separated at about 25% chord downstream of the shock wave and never 
reattached. Distinct periodicity in the computational results was observed and the 
results corresponding to a minimum change in the flow variables between successive 
iteration levels were taken to be the final solution. Figure 5.23 shows the comparison 
between the computational and experimental surface-pressure coefficients. Excellent 
agreement between computational and experimental pressure coefficients is observed 
and the shock location is correctly predicted. The skin-friction distribution on the 
upper surface is portrayed in Figure 5.24 which shows trailing-edge separation. The 
predicted lift, drag and pitching moment coefficients are 0.4395, 0.0176 and -0.0606. 
respectively, which are in close agreement with the corresponding experimental val- 
ues of 0.4378. 0.0163 and -0.0537. The pressure contour plot and the Mach number 
contour plots are shown in Figures 5.25 and 5.26 respectively. 





Figure 5.24: Integrated Technology A 153W: skin-friction coefficient at M 

0.725, Re = 3.0 x 10°, a CO mp = 1-65°. a €X p = 1-65° 
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CHAPTER 6. CONCLUDING REMARKS 


A new nonequilibrium turbulence model has been developed for computation of 
wall bounded two-dimensional turbulent flows. The model, which is patterned after 
the Johnson-King model [12 . is based on the turbulent kinetic energy equation and 
available experimental results for both attached and separated flows. Performance of 
this model was tested by solving the Reynolds time-averaged Navier- Stokes equations 
for various flow conditions around three different airfoil shapes. No changes were made 
in the model parameters while different test cases were being considered. 

Based on the results obtained, both the attached and the shock-induced sepa- 
rated flow cases compare closely with the experimental data. One exception is the 
NACA 0012 separated flow case, where the pressure on the lower surface is slightly 
underpredicted; though not as much as the Johnson-King model [121. It is felt that 
this small difference is due to the modeled form of the turbulence diffusion term. 
Further study on the diffusion term and the inner-layer turbulence velocity scale is 
currently being performed to understand the phenomenon better. Finally, since un- 
steadiness was observed for cases with large regions of separation, the model should 
be tested further by incorporating it into other existing Navier-Stokes codes. 
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APPENDIX A. TURBULENT DIFFUSION 


The turbulent diffusion term is patterned after Johnson and King [12]. As a 
starting point, the bulk convection hypothesis of Townsend [13j, namelv 


dDs 

dy 


Uk y ) 

dy 


( A.l) 


is used. Here V is a lateral convection velocity representing the turbulent transport 
by more energetic large eddies. Along the path of maximum j/t . 


(dDs\ = (y T )m d m (A 2) 

l dy ) m a x dy\y) m 

where the relationship ~t vi/k = a l has been applied. The approximate model for 
the variation of V with y/6 suggested by Johnson and King [12] is utilized here to 
obtain 


L(T\ = ( Y° (A. 3) 

dy \ V / m \0 -tS-ymax/ \ y$n / 

V 0 mentioned in the foregoing equation is the maximum lateral convection velocity 
which is assumed to occur at y — 0-7<5. \q is expected to be proportional to the 


maximum turbulent shear stress, i.e., 
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V o - ^mode^ Tmax ' >1 '^ 

It is to be noted here that the additional restriction on V 0 placed by Johnson and 
King, namely 


i T max J 


dissipation \ 
production } , 


has been dropped here. This is because it is believed here that in the case of equi- 
librium turbulence, the diffusion term is balanced by the convection term in the 
turbulent kinetic energy equation: not that the diffusion term is independently zero. 
Thus equation A. 2 with the aid of A. 3 becomes 


i j ^ 

( dDs \ _ ^ m0( f e / rmr max i Imax 

\ dy J m a i ( 0 . 7<5 - j/mai) ym 
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APPENDIX B. NONEQUILIBRIUM EFFECTS IN THE INNER 

LAYER 


Based on the assumption that the Coles velocity profiles [37; describe turbulent 
mean flows of attached and slightly separated boundary layers, it is possible to deduce 
a relationship between the turbulence velocity scales in the inner and the outer regions 
of such shear layers. Using Hinze's expression for the universal wake function, the 
Coles velocity profiles in defect form is written as 


r 


e - w 


= - si<77l(<7)II ( 1 T COS(7T^)| - In Y 

K L l 0 ) O. 


(B.l 


where 


<7 = 


7 


7 


7 


1/2 


= sign(T w )^~ 

L e 


P w 


Pe 


1/2 


Differentiating B.l with respect to y one obtains, 


du . U 


& 


y I o- i= 

oy 


H) 


+ 


(B.2) 


Since the law of the wall delivers ny | ^ j= crlle — constant , <rUe can be interpreted 
as the turbulence velocity scale VV in the inner layer. The y location for the maximum 
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value of y ^ . as observed from equation B.2, is 0.6466. With these identifications, 
and replacing dujdy with •jJ. equation B.2 may be rewritten as 


*( y ! I )maxp 

1 + 1.82II 


(B.3) 


For equilibrium turbulent flow, (y j \)maxp is proportional to the turbulence 
velocity scale l' 5o in the outer layer. Equation B.3 thus yields a relationship between 
the inner and the outer turbulence velocity scales. 


_ K{constant)V So B 4 . 

s i ~ 1 - 1.8211 

It is now hypothesized that the above relationship holds for nonequilibrium tur- 
bulence as well. Since g e /g can be interpreted as a nonequilibrium correction applied 
to V So . 


K(constant)V Soe g e /g 

V c • — ( D . 0 ) 

S l l-rl.82n 

Finally, writing II in terms of (y j u ; \)max p with the aid of equation B.2 and 
identifying this vorticity function with the equilibrium outer turbulence velocity scale, 
one obtains 

v 'i = V ’i,e Fc (B ' 61 


with 


where Vs- 


f 1 = 


u r ( Pvo / Pe)^^( T w / ' T xv ) [ 1 T wl 1 T w .. 


kFi 


max 


(B.8) 


= crUe has been used. 



